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Mathematical methods together with measurements of single- 
cell dynamics provide unprecedented means to reconstruct in- 
tracellular processes that are only partly or indirectly accessible 
experimentally. To obtain reliable reconstructions the pooling 
of measurements from several cells of a clonal population is 
mandatory. The population's considerable cell-to-cell variabil- 
ity originating from diverse sources poses novel computational 
challenges for process reconstruction. We introduce an exact 
Bayesian inference framework that properly accounts for the 
population heterogeneity but also retains scalability with respect 
to the number of pooled cells. The key ingredient is a stochastic 
process that captures the heterogeneous kinetics of a popula- 
tion. The method allows to infer inaccessible molecular states, 
kinetic parameters, compute Bayes factors and to dissect intrin- 
sic, extrinsic and technical contributions to the variability in the 
data. We also show how additional single-cell readouts such as 
morphological features can be included into the analysis. We 
then reconstruct the expression dynamics of a gene under an 
inducible GALl promoter in yeast from time-lapse microscopy 
data. Based on Bayesian model selection the data yields no 
evidence of a refractory period for this promoter. 

Introduction 

Statistical inference of unobserved molecular states and 
unknown parameters that characterize an intracellular 
process is instrumental for the advance of quantitative bi- 
ology. In particular, single cell assays provide an unprece- 
dented source of data to perform this inference. They 
provide access to the possibly significant stochastic na- 
ture of those processes but also expose the considerable 
heterogeneity in clonal populations of cells. 

From the viewpoint of inference, two classes of single- 
cell data may be distinguished. The first are popu- 
lation snapshot data, provided for instance by cytom- 
etry techniques [H H] |3] or FISH (fluorescence-in-situ- 
hybridization, "5^, where - when measured over time 
- any temporal correlation of a cell's dynamics is nec- 
essarily lost. The second class are time-lapse single-cell 
data f6^, T], where the dynamics of individual cells can 
be followed over several time points. Such a recording 
can be considered a partially observed sample path of the 
underlying process and thus, not only gives us a handle 
on the marginal distribution at each time point but also 
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on how the temporal correlation of the process looks like. 
Evidently, such information is extremely helpful for infer- 
ence and renders live-cell measurements superior in this 
respect - taking aside its limitation in terms of through- 
put compared to flow cytometry, for instance. 

Recently, we proposed a complete inference framework 
for single-cell snapshot data [Ij. In particular, cell-to-cell 
variability was mathematically accounted for and time- 
lapse flow-cytometry data was exemplarily used to infer 
kinetic parameters of a gene expression system in yeast. 
In this work we lay out a corresponding inference scheme 
for time-lapse live-cell data. A very different approach 
needs to be followed here in order to account for the sam- 
ple path nature of the acquired data. 

Several inference techniques for stochastic chemical ki- 
netics based on single-cell time-lapse data have been pro- 
posed recently [8[ [H [TOl [Til HI] • Focus was put on the in- 
ference from the noisy observation of a single sample path, 
because the extension to multiple sample paths is compu- 
tationally straightforward if one assumes no heterogene- 
ity. In practice, pooling of several single-cell recordings is 
necessary to obtain reasonable estimates, which in turn 
generates difficulties due to the extrinsic contribution to 
the observed heterogeneity of a population [131 El fT5] . 
Examples of such extrinsic factors are cell-cycle stage or 
translation efficiency of the cells. The effect of those fac- 
tors on the dynamics of the actual process under study 
has been investigated in fi6\ [TT] . First attempts have 
been made in [181 HI to devise inference techniques that 
account for and estimate such heterogeneity. 

Imaging-based single-cell techniques can additionally 
capture morphological features of cells such as their vol- 
ume or shape. While those were shown to correlate 
well with extrinsic factors [111 [15] , previous inference ap- 
proaches do not offer a principled way to exploit such 
readouts. 

By addressing extrinsic variability, the present work 
goes significantly beyond previous inference efforts. 
Based on the innovation theorem for counting processes 
[2QJ , we obtain a marginal stochastic process by integrat- 
ing the heterogeneous Markov process over the extrinsic 
factors as well as fixed kinetic parameters with prior un- 
certainty. This results in a scalable inference scheme, the 
parameter dimension of which is independent of the popu- 
lation size. In conjunction with the time-lapse microscopy 
measurements, the method allows us to reconstruct the 
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dynamics of an inducible gene expression system in yeast 
from pooled time-lapse microscopy data. 

Results 

Mixed- effect Models and Measurements. We con- 
sider a continous-time Markov chain (CTMC) describing 
the time-evolution of a well-mixed reaction system with 
d chemical species and u reaction channels. The random 
state vector X{t) collects the copy numbers of species at 
time t and random paths thereof on the time interval [0, t] 
are denoted X. Subsequently, we will use the convention 
to denote a random quantity by an upper case letter and 
its realization by the corresponding lower case letter. The 
propensity functions of the reactions are assumed to be 
of the form hj{x) = CjQjix) with Cj the stochastic rate 
constant and gj a function of the current state x. More- 
over, we denote the set of all stochastic rate constants as 
C = (Ci, . . . , Cu)' Following a Bayesian approach, the 
random parameters C get equipped with a prior believe 
through a probability distribution p{c). Note that besides 
the reaction rate constants, a single-cell stochastic model 
of a cellular process may additionally be parametrized 
in terms of its initial conditions, size of the reaction com- 
partment, the state of the molecular environment in which 
the process is embedded and so forth. We denote such a 
complete parametrization of the Markov chain as B. 

Apart from intrinsic fluctuations of the stochastic pro- 
cess, extrinsic variations will render process dynamics dif- 
ferent across cells in a population even in the absence of 
intrinsic fluctuations - suggesting the use of individual 
Markov chains for each cell. Hence, we allow some of the 
process parameters to vary across cells. We denote those 
as extrinsic factors Z C Q. They can either be assumed 
constant (TJ [5] or dynamically changing [T71 [16] . Extrin- 
sic factors include rate constants of aggregated multi-step 
reactions that dependent on the presence of intermediates 
(e.g. translation rate), initial conditions or similarly total 
copy numbers of conserved proteins, ribosome numbers, 
cell volume or cell-cycle stage - to name just a few. Here, 
we take a macroscopic view and address the evident het- 
erogeneity in a cell population that is believed to roughly 
stay unchanged during the course of an experiment. In 
turn, we neglect the microscopic picture [IX J6j that some 
extrinsic factors, on top of being different across cells, also 
fluctuate temporally and modulate the considered pro- 
cess. For instance, ribosomes are themselves subject to 
intrinsic fluctuations. Hence, we assume that those fluc- 
tuations are small compared to the constant offsets across 
cells and to the intrinsic fluctuations of the actual process 
under study. On the other hand, some extrinsic factors, 
such as the cell-cycle stage may not fluctuate but rather 
vary deterministically and slowly with respect to the ex- 
perimental time span [HI [HJ [53j , such that quasi-steady 
state arguments can support a time-invariance assump- 
tion. 

In statistics the resulting heterogeneous population 
models are known as mixed- effect models, where the com- 
plement subset S = Q \ Z are shared among cells and 
correspond to the fixed effects, whereas the extrinsic fac- 
tors refer to the random effects. Examples in S are rate 
constants of elementary reactions, such as protein associ- 



ations that are determined by the biophysics of the two 
proteins and should not vary across cells. Following the 
time-invariance of extrinsic factors Z, we assign to them 
a probability distribution p{z \ a), where the extrinsic 
statistics a is a set of parameters specifying the shape 
of this distribution and hence the population's hetero- 
geneity. Given a population of M cells, the m-th cell's 
state is then described by a conditional Markov chain 

Xm I ^Z'^^S). 

Experimentally we can retrieve a corrupted abundance 
measure for a small subset of molecular species for M 
different cells at each measurement time ti. The asso- 
ciated acquisition error is characterized by a conditional 
measurement density p{yl^ \ with x]^ = x^{ti) 

and uj an unknown distribution parameter such as the 
acquisition noise variance. Furthermore, we define the 
state trajectory of cell m between the l-th and the k-th 
measurement time as X]!]. and denote by the cor- 
responding set of measurement^. Morphological features 
are accounted for by introducing morphological covariates 
and hypothesizing a statistical dependency between 
these covariates and the extrinsic factors Z^. This depen- 
dency is described by a conditional density p{v^ \ z^^b)^ 
with b a set of shape parameters characterizing this condi- 
tional density. A schematic representation of the resulting 
mixed-effect Markov model is given in Figure la and the 
corresponding hierarchical Bayesian network is depicted 
in Figure lb. 

Heterogeneous Kinetics. According to the above 
mixed-effect model, the dimension of the parameter space, 
on which we aim to perform inference subsequently, is 
scaling with the number of considered cells M. Hence, a 
marginalized process model, where the unknown extrin- 
sic factors are integrated out, appears critical to achieve 
practical tractability. Formally, the dynamic state of each 
cell m is then described by one and the same marginal 
process | {S,A) := X \ {S,A) with path density 

p(x I s, a) = J p(x I s)p{z I a)dz. For the sake of 
illustration, let us assume & = C and Z to be one- 
dimensional, namely the stochastic rate constant of reac- 
tion j that hence randomly varies across the population. 
Based on the innovation theorem for counting processes 
|2U| , we show that the marginalized dynamics again fol- 
lows a jump process, where the propensity corresponding 
to the extrinsic parameter can be generally written as 



/i,(x,t)=E[Z|x,a]^,(x(t)), 



(1) 



with E [Z I X, a] as the conditional expectation of Z given 
a complete sample path x and the extrinsic statistics a. 
For instance, if Z follows a Gamma distribution p{z \ 
a) = with a = we have that p{z | x, a) = 

g{a^rj,(3^ J^gj{x{r))dr) [Ml US] and hence. 



hj{yi,t) = 



a + ri 



-9j{<t)). (2) 



le number of occur 
j in X. For CTMCs, ([1]) can be rewritten [26j in terms 



where rj denotes the number of occurrences of reaction 



■•^Note that in contrast to Yj^, the state trajectory XJT^^ denotes 
a continuous-time sample path between ti and ■ 
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Extrinsic factors - latent 





Figure 1: Modeling heterogeneous microscopy data, (a) Schematic generative model of the experimental data. On top of intrinsic 
fluctuations, extrinsic factors and their morphological covariates render individual cells different. Recorded intensitiy levels are hence 
characterized by a parameter set S that is shared across cells and a set of individual (i.e. extrinsic) parameters Z*. (b) Corresponding 
Bayesian mixed-effect model. Nodes denote random variables and statistical dependency (and causality) is indicated by directed 
edges. Gray-shaded nodes correspond to experimentally accessible quantities, whereas all other nodes refer to latent variables. 
Extrinsic factors are drawn from a common distribution parametrized by A. Note that for clarity, the nuisance parameters Q and 
B defined in the text are not shown. 



of the moment-generating function of Z \ A; mathemat- 
ical proofs and derivations for the general case of multi- 
variate extrinsic variables are provided in section S.l in 
the Supporting Information. If morphological covariates 
are included into the analysis, the marginalization can 
be performed analogously by constructing a marginal- 
conditional process \ (6', A, 5,V^) with path den- 
sity j9(x^ I s^a^b^v^) (see section S.l in the Supporting 
Information). With the above marginalization, the ran- 
domness over extrinsic factors is exposed by the stochas- 
tic process itself which now shows a direct dependency 
on the extrinsic statistics A. Note that as a consequence, 
the Markov property is lost and moreover, the propen- 
sity hj{'K^t) becomes explicitly time- varying. However, 
when augmenting the state space by the summary statis- 
tics T(x) = ( Tj, Jq gj{x{T))dr ) of the path x, the Markov 
property is recovered and hence, stochastic simulation can 
efficiently be performed using available methods [27], see 
section S.1.1 in the Supporting Information). On top 
of several computational advantages, the marginal pro- 
cess provides an important means to properly describe 
stochastic biomolecular dynamics in presence of extrin- 
sic variability. More specifically, it reflects the dynam- 
ics of a cell with unknown extrinsic factors and hence. 



coincides with typical real- world scenarios. In contrast 
to traditional stochastic models that rely on full para- 
metric knowledge of the considered process, the marginal 
jump process can take into account uncertainty over ex- 
trinsic factors and thus, provides a principled framework 
to study kinetics in a heterogeneous population. 

Statistical Inference. Before we outline the inference 
procedure in more detail, we discuss an important impli- 
cation of the above marginalization approach. In particu- 
lar, it follows that any kind of parametric uncertainty, lin- 
early entering the propensities can be integrated out and 
directly encoded into the stochastic process. This seems 
particularly useful from a Bayesian viewpoint, where pa- 
rameters are characterized by prior uncertainty. For in- 
stance, if the i-th kinetic parameter Ci is associated with 
a known prior distribution, e.g. p{ci) = Q{ai,f3i), the 
marginal propensities are obtained in analogy to (|2]). 
Marginalization with respect to every kinetic parameter 
in conjunction with a Bayesian filtering approach leads to 
the following inference scheme, which we refer to as dy- 
namic prior propagation (DPP). The goal of the scheme 
is to compute the marginal posterior distribution 
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where {y^N^v^) denotes the tuple of measurements sample from (|3]). However, by the curse of dimension- 
available for the m-th cell. In principle, a Markov chain ality such approaches are prohibitive in terms of accep- 
Monte Carlo (MCMC) scheme can be applied directly to tance ratios when considering multi-dimensional dynam- 
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ics over multiple measurement time points. According 
to a Bayesian filtering approach, the posterior distribu- 
tion is determined by applying Bayes' law recursively over 
time. Consequently, the original sampling problem breaks 



up into a sequence of sub-problems with significantly re- 
duced dimensionality. The posterior distribution at time 
ti is then given by 



p (xl, ,...,^^„a,b,u;\ {yl,, v'), . . . , {y^^,v^)) ^ 



M 



n p{yr I xr,io)p{,dr-i:i i xr.„a,b,v"',T{^^i_,)) 



(4) 
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Although analytically intractable, we can sample from 
(HI) using a Metropolis-Hastings sampler, for instance. Re- 
cursive sampling approaches inherently suffer from sam- 
ple degeneracy as soon as constant parameters are esti- 
mated in addition to the dynamic states. While a ma- 
jority of the parameters is integrated out thanks to the 
marginal process framework, A, B and Q remain in the 
model. A standard approach to avoid such degeneracies 
is to apply an invariant kernel to the static parameters at 
each time instance, such as to maintain diversity among 
parameter samples [28j. Here, we use a Metropolis- 
within-Gibbs MCMC scheme, where the latent space is 
further divided into blocks which are successively resam- 
pled. This requires sampling from the full conditional dis- 
tributions, which can be determined using the notion of 
Markov blankets f535 . A full description of the algorithm 
can be found in section S.2 in the Supporting Information. 
Note that posterior distributions over Z^, . . . , and S 
are not directly computed by the marginalized inference 
scheme. However, they can be easily reconstructed via 
the law of conditional probability, such as described in 
section S.2. 4 in the Supporting Information. 

Synthetic Data: Two-State Gene Expression 

Model. We first studied the proposed inference frame- 
work using synthetic data of a simple two-state gene ex- 
pression model [4J given in Figure 2a under realistic mea- 
surement conditions. We assume that the target gene 
can be activated upon application of an exogeneous sig- 
nal - for instance via stress-induced translocation of cer- 
tain transcription factors. Extrinsic variability was simu- 
lated by introducing a Gamma-distributed variability in 
the translation rate C5. Data was collected for M = 20 
cells, on which we applied DPP using 10000 samples per 
time instance (for further details, see caption of Figure 
2a and Figure 2b). The inferred posterior distributions of 
the kinetic parameters, the extrinsic statistics as well as 
the acquisition noise parameter are depicted in Figure 2b 
(for a density plot of the morphological shape parameters 
see Figure S.5 in the Supporting Information). Further- 
more, state inference with respect to mRNA and protein 
levels is demonstrated in Figure 2c using two exemplary 
cells with different translation efficiencies. 

Importantly, the inferred posterior distribution over the 
latent states can be used to reconstruct switching events 
of the target gene. We remark that a simple mode detec- 
tion on the marginal distribution (time-point-wise) does 
not yield a valid sample path and hence cannot be used to 
extract timing statistics, such as gene on/off times. Con- 
sequently, we determined the most likely posterior (MAP) 
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Figure 2: Parameter and state inference using synthetic data, 
(a) Two-state gene expression model. The gene activation event 
is controlled by a time-varying rate u{t) and subsumes transcrip- 
tion initiation events such as recruitment of transcription factors, 
RNA polymerase and possibly chromatin remodeling. Upon ac- 
tivation of the gene, mRNA can be transcribed to further yield 
new proteins. All reactions are modeled according to mass-action. 
Altogether, the model comprises four species and six reactions, 
whereas we assume a Gamma-distributed heterogeneity in the 
translation efficiency (i.e., C5) with distribution parameters a and 
p. We assume that = 10 noisy (log-normally distributed with 
unknown scaling parameter) measurements of the protein abun- 
dance can be obtained at equally spaced time points within a 
total interval of roughly 40 minutes, (b) Parameter inference 
from protein time series using M = 20 cells. Inference results 
are shown for the three kinetic parameters (c2, C3 and C4), the 
extrinsic statistics a and /3 (describing the heterogeneity over C5) 
and the scaling parameter cj of the acquisition noise. The se- 
quential MCMC scheme was performed using 10000 samples per 
time instance, (c) Inferred mRNA and protein abundance. The 
5%- and 95%-quantiles of the inferred state distributions as well 
as their mean values, i.e., minimum mean square error (MMSE) 
estimates are shown for two representative cells with different 
translation efficiencies; the gray-shaded area indicates the win- 
dow of induction. The thick violet line illustrates the assumed 
Gamma distribution over the extrinsic factor, i.e., the translation 
rate. 
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Figure 3: Reconstructed gene expression dynamics for a 
double-pulsed induction from noisy protein measurements (cir- 
cles). Shaded areas denote 5%- and 95%-quantiles of the pos- 
terior distributions over states and the true values (solid), the 
MMSE estimates (dashed) and the maximum a posteriori (MAP) 
paths (blue) are shown. Also shown are the posterior probabilities 
for the gene to be inactive or active, Po{t) and Pi{t), respec- 
tively; black and white coloring corresponds to probability one 
and zero, respectively. Furthermore, we computed the marginal 
MAP (MMAP) state estimate (green) based on Po{t) and Pi(t). 

switching path of the gene from which timing statistics 
may be recovered. This approach to reconstruct switching 
events is exact and contrasts other approximate schemes 
[6], the vahdity of which was not demonstrated on syn- 
thetic data where the actual switching path is known. 
In general, the inverse problem of reconstructing a gene 
switching path from the slow protein dynamics is consid- 
erably ill-posed and we expect the posterior density over 
the switching paths to be close to degenerate. However, 
for the simulation study. Figure 3 indicates that accurate 
detection of the gene state is indeed possible within the 
realistic scenario considered here. For this test we simu- 
lated a pulsed induction of gene expression. 

Experimental Data: /3-Estradiol-Induced Gene 
Expression. The DPP-based inference was used to re- 
construct the expression dynamics of an artificially con- 
trolled gene expression system in yeast. A widely used 
system to control the expression of genes under a GALl 
promoter in Saccharomyces cerevisiae is based on the 
hormone-dependent activation of the chimeric transcrip- 
tion factor GAL4DBD.ER.VP16 (GEV) [3D1 |3l]. GEV 
consists of a strong transcriptional activator, made by 
fusing the GAL4 DNA binding domain (GAL4DBD) with 
the hormone-binding domain of the human estrogen re- 
ceptor (ER) [30j and the transcription activating domain 
of the herpes simplex virus protein VP 16 [32]. In its in- 
active state, GEV associates with the Hsp90 chaperone 
complex and resides in the cytoplasm. Upon addition of 
the exogenous hormone /3-estradiol to the extracellular 
medium, /3-estradiol diffuses through the cell membrane 
and binds to the GEV's ER. Thereby, Hsp90 disassociates 
from the complex and active GEV translocates to the nu- 
cleus where it's GAL4DBD recognizes and binds to GAL 
promoter regions. VP16 then activates transcription of 
the downstream gene. 

We engineered a strain that allows a combined readout 
of GEV translocation and /3-estradiol-induced gene ex- 



pression. A fully functional GEV-mCherry construct in 
combination with a nuclear marker allows to compute a 
ratio of nuclear to cytoplasmic GEV, serving as our model 
input. In the same strain, a destabilized [33j version of 
the Venus fluorescent protein (Y- Venus) was placed un- 
der control of a GALl promoter, which allows a more 
accurate tracking of the gene expression dynamics. For 
a detailed description of the strain see S.4.1 in the Sup- 
porting Information. 

The fluorescence microscopy experiments were per- 
formed using a flow chamber that allowed us to rapidly 
exchange the extracellular media and apply a 30min pulse 
of 50nM /3-estradiol. The time-lapse microscopy movies 
were automatically analyzed [3lj to quantify the change 
in nuclear localization of the GEV-mCherry protein and 
the fluorescent levels of the Y- Venus protein in individ- 
ual cells (see S.4.5 in the Supporting Information). A 
calibration curve was used to map recorded intensities to 
total protein abundances (see S.4.3 in the Supporting In- 
formation). Based on a one-sided Kolmogorov-Smirnov 
test we determined and corrected for a time delay in the 
protein measurements, that arises from an aggregation 
of unmodeled sequential events, such as mRNA export, 
post-transcriptional/translational modifications and re- 
porter maturation (see S.4.3 in the Supporting Informa- 
tion). We then used M = 20 single cell trajectories of 
Y- Venus abundance and the average GEV translocation 
within the subsequent analyses. 

Modeling GALl Y- Venus Expression and Acquisi- 
tion. We determined empirical evidence for three mod- 
els of eukaryotic gene expression by Bayesian model selec- 
tion. Next to the canonical two-state model [l] (see Figure 
4a) we considered a model with a third refractory state 
[71 [6] (see figure S.6b in the Supporting Information) and 
a three-state variant of a model proposed in [35], where 
the initiation-complex assembly is followed by a slow acti- 
vation step representing either RNA polymerase (RNAP) 
binding or chromatin remodeling (see figure S.6c in the 
Supporting Information). Bayes factors of roughly 2 and 
43 dB, respectively, in favor of the two-state model were 
found. Bayes factors were also computed for two com- 
peting measurement noise models (i.e., i.i.d. normal and 
log- normal). Strong evidence was found for log- normally 
distributed measurement noise, which was reflected by a 
highly decisive Bayes factor (> lOOdB). 

GALl Transcriptional Dynamics Exhibit Mild 
Bursting. The posterior distributions over the kinetic 
parameters were computed via DPP using 20000 samples 
during the first - and 10000 samples during the subsequent 
time iterations (see Figure 4b). State reconstruction and 
gene activity detection was performed as described pre- 
viously. Figure 5a shows the reconstructed dynamics for 
two cells with different Y- Venus abundances. The order 
of the estimated mRNA half-life of around 10 min and 
synthesis rate of 6 molec/min are inline with previous 
findings [36j. The specific value of the latter is above 
most reported rates for constitutively expressing genes 
[36] which appears consistent with the fact that inducible 
GAL-promoter driven genes show high expression levels 
that can even lead to toxicity and protein-aggregation 
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Figure 4: Parameter inference for a model of /3-estradiol- 
induced Y-Venus expression, (a) Stochastic model of /3-estradiol- 
induced gene expression in yeast with extrinsic noise introduced 
by a variable translation efficiency C5 (blue highlighted), (b) Pos- 
terior distributions over unknown parameters. Based on 20000 
samples for the first time-iteration and 10000 samples for the 
subsequent iterations dynamic prior propagation was performed 
with respect to all kinetic parameters (ci, . . . , C4, ce), the acqui- 
sition noise parameter cj as well as the extrinsic statistics {a, 13) 
and the morphological shape parameters (yo,0) that characterize 
the translation rate C5. The marginal posterior for the gene-on 
rate is shown for au, with u as the temporal average over the 
modulating GEV intensity. 



when constitutively induced ^37J. This synthesis rate to- 
gether with the length of 850 bp for the Y-Venus protein 
and a reported elongation speed of 2 kb/min [38] for GAL- 
driven genes indicates that there need to be at least four 
RNAPs on average on the gene. The timing statistics 
obtained from reconstructed gene-switching paths yields 
that for successful initiations on average around 2.5 tran- 
scripts per active gene state are produced, suggesting that 
transcription re-initiation and thus mild bursting takes 
place in this expression system. 

Noise Contributions in GALl Y-Venus Expres- 
sion. Figure 5a also indicates to which extent a cell's 
expression level is explained by extrinsic and intrinsic fac- 
tors. Although cell 1 shows mRNA levels similar to cell 
2, the former expresses significantly more Y-Venus due to 
a larger translation rate C5. By forward simulating the 
inferred model we dissected and quantified the different 
sources of variability in the measured Y-Venus abundance 
using the law of total variance (see Methods). In partic- 
ular, we separated intrinsic, extrinsic and technical con- 
tributions to the overall variability, which are shown in 
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Figure 5: Reconstruction of heterogeneous Y-Venus dynam- 
ics, (a) Inferred dynamics for two exemplary cells showing differ- 
ent Y-Venus abundance. The temporal GEV-induction is shown 
as an intensity map, whereas white and black coloring denote 
minimal and maximal abundance, respectively. The inlet scat- 
ter plots indicate the inferred expression regime, where each dot 
represents a single cell and the arrows point towards the corre- 
sponding cell; the x- and y-axes correspond to the logarithm of 
the temporal mean of Y-Venus and the MMSE estimate of C5, 
respectively. The panel below shows the inferred mRNA dynam- 
ics. Therein, the inlet with the transcripts per on cycle distri- 
bution over all posterior paths shows a mean around one (gray 
dashed) and mean around 2.5 taking successful initiation events 
(red dashed), (b) Sources of cell-to-cell variability in Y-Venus 
expression. The inferred model was used to simulate M = 10000 
trajectories which were used to compute the squared coefficient of 
variation (SCV) of the Y-Venus abundance. The total SCV was 
decomposed into technical, intrinsic and extrinsic components 
(see Methods), (c) Dependency between volume increase and Y- 
Venus abundance. Multiple cell's intensity trajectories, their vol- 
ume increase and translation efficiency (i.e., C5) were computed 
via forward-simulation of the inferred model. The plot shows the 
Y-Venus abundance at time T = 200min versus the volume in- 
crease from T = Omin until T = 200min for the predicted (blue) 
and experimental (red) data. The statistical dependency between 
the volume increase and the translation efficiency is indicated by 
the dashed iso-lines (gray). 



Figure 5b. In the technical contribution any systematic 
bias introduced, for instance by the image segmentation 
algorithm, is not considered. The inferred model predicts 
that the variability in Y-Venus expression is dominantly 
driven by extrinsic factors. Although the technical con- 
tributions appear small compared to the total population 
variability, each intensity trajectory is characterized by 
significant measurement uncertainty (such as indicated 
in Figure 5a). 

Morphological Features and Extrinsic Variability. 

On top of the protein measurements, the inference pro- 
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cedure allows to incorporate additional single-cell read- 
outs such as morphological features. More specifically, it 
is able to quantify statistical dependencies between such 
readouts and a population's extrinsic factors, offering an 
efficient way to generate and validate biological predic- 
tions. Here we hypothesized a dependency between vol- 
ume increase during the observation time interval and 
the extrinsic factor (i.e., the translation efficiency) and 
quantified it using the proposed inference algorithm. The 
found covariation is depicted in Figure 5c. Consistent 
with [H] we find that volume increase correlates with 
translation efficiency (r = 0.56, p = 0.01) but that it does 
not explain all extrinsic variability present in the inten- 
sity trajectories. Hence, this provides more evidence for 
the fact that simple normalization through morphological 
features (e.g. forward scattering in ffow cytometry data 
[1]) can not sufficiently correct for extrinsic variability in 
data. 

Discussion 

Inference of stochastic dynamics of a cellular process from 
live-cell recordings of a single cell is inherently ill-posed. 
In contrast, pooling together of even a few cell record- 
ings (e.g. around 10-20) significantly improves the accu- 
racy of the inference and can even resolve practical non- 
identifiabilities (see S.3.4 in the Supporting Information). 
However, the large degree of extrinsic variability being 
evident in such data induces another layer of complexity 
where straightforward approaches suffer from the curse of 
dimensionality. 

The proposed framework rests upon a recursive infer- 
ence scheme, whose strength is achieved by integrating 
out extrinsic factors and uncertain kinetic parameters 
from the process dynamics. First, kinetic parameters are 
no longer sampled along with the dynamic states result- 
ing in a variance reduction of the desired posterior statis- 
tics (see S.3.1 in the Supporting Information). Second, 
in the context of cell-to-cell variability, the marginaliza- 
tion yields scalability with respect to the cell population 
size and therefore appears advantageous to previous ap- 
proaches that require intermediate sampling of extrinsic 
factors [18j. The approach is conceptually related to the 
rao-blackwellized particle filter [39j, where a set of lin- 
ear dynamic states is marginalized out in order to reduce 
the dimensionality. While those approaches are limited to 
a special class of system dynamics, DPP applies to any 
non-linear reaction network with linearly-parametrized 
propensities, such as obtained by mass-action principles. 
Apart from inference, the marginalized process is impor- 
tant in itself to properly describe the dynamics of hetero- 
geneous cell populations. We remark that the presented 
inference framework is exact in the sense that for the limit 
of infinite samples the computed distribution converges to 
the correct posterior smoothing distribution. 

The simple two-state gene expression model from Fig- 
ure 2a demonstrates correctness of the novel inference 
framework. Although only a few measurements of the 
protein abundance were available, all parameters (see Fig- 
ure 2b) as well as gene-switching sequences (see Figure 3) 
were accurately inferred from the data. 

When applied to the engineered /3-estradiol-induced 



gene expression system in yeast, the inference scheme to- 
gether with performed protein calibration measurements 
offers estimates in terms of absolute quantitation for ki- 
netic parameters, unobserved molecular states and popu- 
lation heterogeneity. Hence, it allows to estimate several 
quantities and their uncertainty at once without resorting 
to dedicated experimental techniques H] . The results 
are given for the two-state gene expression model because 
Bayes factors indicated weak empirical evidence for more 
complex models. In particular, no refractory gene state 
was evidenced by the data. 

By pooling heterogeneous single-cell traces we can 
model-based dissect the different contributions to cell-to- 
cell variability - something that traditionally requires ex- 
periments such as two-color assays fi3\ HP] . Inline with 
previous studies on GAL-driven genes [35l we find 
that extrinsic noise is the dominating source of cell-to- 
cell variability for such genes. Moreover, in the course of 
induction we observe, consistent with Poissonian noise, 
a characteristic decrease of the intrinsic noise compo- 
nents accompanying the increase in mean expression level. 
The magnitude of technical or measurement noise was 
estimated to be in the same order than that of intrin- 
sic noise. Bayesian model selection for the measurement 
noise model provided evidence in favor of log-normally 
distributed noise with respect to the normal distribution. 
This indicates that measurement errors scale with the 
measured fiuorescent intensity, such as incurred in image 
segmentation. Morphological features [19j in addition to 
intensity trajectories can significantly increase the pre- 
dictive power of computational models with respect to 
extrinsic factors. The probabilistic dependency between 
translation efficiency and volume increase extracted by 
the algorithm is coherent with earlier findings [jA^ that 
both quantities are positively correlated. 
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Methods 



References 



Supporting Information. A PDF document with 
supplementary derivations and details is available at 
http:/ /www. bison. etliz.cli/bison/SI_Pooling/, 



Experimental Procedures. Measurements were per- 
formed in an incubation chamber at 30°c using an au- 
tomated epi-fluorescence microscope (Eclipse Ti, Nikon 
Instruments) with a 60x oil objective and specific 
(CFP/YFP/mCherry) excitation and emission filters. 
Images were acquired at multiple positions using a motor- 
ized XY- stage, while the focal plane was maintained using 
a hardware autofocus system. Cell chambers (/i-siide vi°-^, 
ibidi) were treated with filtered solution of Concanavalin 
A dissolved in PBS (img/mi) for somin and subsequently 
rinsed with PBS. Single colonies of the respective yeast 
strain were picked and inoculated in synthetic (SD) 
medium. The overnight saturated cultures were then di- 
luted and grown in log-phase for at least two doubling 
times (>4h), again diluted (OD600 0.01), briefly sonicated 
and loaded into the cell chambers. Pulse experiments 
were performed by switching between two hydrostatic 
pressure driven liquid flows (see section S.4 in the Sup- 
porting Information). A calibration curve was obtained 
by measuring fluorescence intensities for several proteins 
of known average abundance, tagged with a Venus fluo- 
rescent reporter. For a graphical illustration and further 
details see S.4. 3 in the Supporting Information. 



Inference Algorithm. A full technical description of 
the sampling-based inference algorithm can be found in 
section S.2 in the Supporting Information. Empirical 
model evidences - and Bayes factors are directly calcu- 
lated by the algorithm and do not require further compu- 
tations. The corresponding formulas and their derivations 
are given in S.2. 6 in the Supporting Information. A cor- 
responding Matlab toolbox with a simple user interface is 
made available under the GNU public license. 



Model-Based Noise Decomposition. The law of to- 
tal variance is applied to dissect the total variability into 
intrinsic, extrinsic and technical contributions. In partic- 
ular, it holds that 



SCV[Yi] = 



E [E [Var [Yi \ Xi] \ Z]] 



technical 

E[Var[E[l; | Z]] 



Var [E [E [Yi \ Xi] \ Z]] 



intrinsic 



For the decomposition, the model parameters were set to 
their inferred MMSE estimates and the individual quanti- 
ties were obtained via forward simulation. Further details 
can be found in section S.2. 7 in the Supporting Informa- 
tion. 
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